Strong-coupling perturbation theory for the extended Bose-Hubbard model 
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We develop a strong-coupling perturbation theory for the extended Bose-Hubbard model with on-site and 
nearest-neighbor boson-boson repulsions on {d > l)-dimensional hypercubic lattices. Analytical expressions 
for the ground-state phase boundaries between the incompressible (Mott or charge-density-wave insulators) and 
the compressible (superfluid or supersolid) phases are derived up to third order in the hopping t. We also briefly 
discuss possible implications of our results in the context of ultracold dipolar Bose gases with dipole-dipole 
interactions loaded into optical lattices. 
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I. INTRODUCTION 

Ultracold atomic physics in optical lattices has created a 
new experimental arena where many simple model Hamiltoni- 
ans can be constructed and 'simulated' experimentally [IJ]. To 
date, the most successful efforts have been with bosonic atoms 
on optical lattices 10, H, 0,111]. Here, when the single-particle 
bands of the optical lattice are well separated in energy, the 
boson-boson interaction is much smaller than that separation, 
and the particle filling is not too high, the system is described 
well by the single-band Bose-Hubbard (BH) model. This 
model is the bosonic generalization of the Hubbard model and 
was introduced originally to describe "^He in porous media or 
disordered granular superconductors |6]. The superfluid phase 
of bosonic systems is well described by weak-coupling theo- 
ries, but the insulating phase, where there is a gap to particle 
excitations with a uniform (integer) filling of the bosons on 
each lattice site, is a strong-coupling phenomenon that only 
appears when the system is on a lattice. This Mott insulator 
phase is incompressible and hence occupies a finite area in the 
parameter space of the chemical potential and the hopping. It 
has a transition from the incompressible phase to a compress- 
ible superfluid as the hopping or chemical potential are varied. 
The on-site BH model has been studied extensively, and the 
strong-coupling perturbation theory approach has been shown 
to be quite accurate in determining this phase diagram of the 
system 13 iSfl. 

Recently, experimental progress has been made in con- 
structing ultracold dipolar gases of molecules, namely K-Rb 
molecules, from a mixture of fermionic ""^K and bosonic '^^Rb 
atoms 1 9, 10]. In this case, the molecules are fermionic, but 
similar principles will allow one to also create bosonic dipo- 
lar molecules. Future experiments are likely to load these 
bosonic molecules into optical lattices. These systems will 
have a long-range boson-boson interaction mediated by their 
dipole moment, which can be approximated, in some circum- 
stances, by an on-site and a nearest-neighbor repulsion (gener- 
ically, dipole-dipole interactions will be longer ranged than 
just nearest neighbors and also can have directionality due to 
the orientations of the dipoles). The case of an extended BH 
model, where the boson-boson interaction is longer ranged, 
has also been widely studied lUll [H [H [H Si- Inclu- 



sion of a nearest-neighbor repulsion can lead to the formation 
of a charge-density-wave (CDW) phase, where, at half-filling 
for example, one would have a checkerboard arrangement of 
the density in an ordered pattern. This phase is incompress- 
ible with a finite gap to excitations; it also breaks the origi- 
nal translational symmetry of the lattice, forming a new crys- 
talline phase. The CDW phase has generated significant inter- 
est, because it often can become a supersolid prior to becom- 
ing a superfluid as the interactions are reduced. A supersolid 
phase is a (compressible) superfluid that continues to have a 
density modulation (or CDW) present 1 16.1; that is, the su- 
perfluid and crystalline orders co-exist. Interest in supersolid 
physics has increased dramatically since the recent observa- 
tion of supersolid-like behavior in low-temperature He exper- 
iments 1 17]. There is some numerical and theoretical evidence 
that the supersolid phase exists only in dimensions higher than 
one iQllil. 



In this work, we examine the extended BH model with 
on-site and nearest-neighbor boson-boson interactions via a 
strong-coupling perturbation theory in the hopping, plus a 
scaling analysis, which allows us to accurately predict the crit- 
ical point, and the shape of the insulating lobes in the plane of 
the chemical potential and the hopping. We carry the analysis 
out to third-order in the hopping, and we perform the scaling 
theory using the known critical behavior at the tip of the in- 
sulating lobes [which corresponds to the (d + l)-dimensional 
XY model, and is identical for the Mott and CDW phases]. 



The remainder of the manuscript is organized as follows. 
After introducing the model Hamiltonian in Sec. [Ill we de- 
velop the strong-coupling perturbation theory in the kinetic- 
energy term in Sec.|IIIl where we derive analytical expressions 
for the phase boundaries between the incompressible (Mott or 
CDW insulators) and compressible (superfluid or supersolid) 
phases. There we also propose a chemical-potential extrap- 
olation technique based on scaling theory to extrapolate our 
third-order power series expansion into a functional form that 
is appropriate for the Mott or CDW lobes, and compare these 
results with the mean-field ones in Sec. |IV] A brief summary 
of our conclusions is presented in SeclVl 
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II. EXTENDED BOSE-HUBBARD MODEL 

We consider the following extended BH Hamiltonian with 
on-site and nearest-neighbor boson-boson repulsions 

+ ^Vijninj - ^i^n^, (1) 

i,j i 

where <y is the tunneling (or hopping) matrix between sites 
i and j, b\ (bi) is the boson creation (annihilation) operator 
at site i, rii = b\bi is the boson number operator, J7 > is 
the strength of the on-site and Vij is the longer-ranged boson- 
boson repulsion between bosons at sites i and j, and /i is the 
chemical potential. In this manuscript, we assume tij is a 
real symmetric matrix with elements ty — t for i and j near- 
est neighbors and otherwise and similarly for Vij (equal to 
V > for i and j nearest neighbors and zero otherwise), and 
consider a (d > l)-dimensional hypercubic lattice with M 
sites. Note that we work on a periodic lattice with no external 
trap potential. 

We also assume U > zV where z = 2dis the lattice coor- 
dination number (number of nearest neighbors). In this case, 
the boson occupancy of the nearest-neighbor sites in the CDW 
phase can only differ by one. For instance, the first CDW 
phase is such that every other site is occupied by one boson 
and the remaining sites are left unoccupied. When U < zV, 
additional CDW phases can be present in the phase diagram. 
For instance, a CDW phase in which every other site is oc- 
cupied by two bosons and the remaining sites are left unoc- 
cupied is energetically more favorable than a Mott phase in 
which every lattice site is occupied by one boson. Our results, 
with minor changes, can also be used to analyze these addi- 
tional CDW phases if desired, but more work would be needed 
to examine other types of CDW order, like columnar (stripes) 
and so on, which can arise from longer-range interactions. 



from uq to no + 1 when /i = Una + 0^. For instance, the 
ground state is a vacuum with riQ = for /i < 0; it is a Mott 
insulator with uq — 1 for < /i < ?7; it is a Mott insulator 
with Uq — 2 for U < fj, < 2U, and so on. 

When V 0, the ground state has an additional CDW 
phase which has crystalline order in the form of staggered 
boson densities, i.e. (rii) — Ua and (rij) = rifc for i and j 
nearest neighbors. Therefore, to describe the CDW phases, it 
is convenient to split the entire lattice into two sublattices A 
and B such that the nearest-neighbor sites belong to a different 
sublattice (a lattice for which this can be done is called a bi- 
partite lattice — we assume the number of lattice sites in each 
sublattice is the same here). We assume that the boson occu- 
pancies of the sublattices A and B are ria and rib, respectively, 
such that Ua > rn,. We remark that the Ua = = no states 
correspond to the Mott phase. It turns out that the ground- 
state energy of the {ria = no + 1, rib = no) state is degenerate 
with those of the (ua = no, nf, = no) and (n^ = no + 1, n?, = 
?io + l) states at /i = Uno + zVno and = Uno + zV{no + l), 
respectively. This means that the chemical potential width of 
all Mott and CDW lobes are U and zV, respectively, and that 
the ground state alternates between the CDW and Mott phases 
as a function of increasing /i. For instance, the ground state 
is a vacuum {ua = 0, n^ = 0) for /i < 0; it is a CDW with 
{ua = 1, nf, = 0) for < /i < zV; it is a Mott insulator with 
{ua = l,nb = 1) for zV < < U + zV; it is a CDW with 
{ria = 2,nb 1) for U + zV < (U < ?/ + 22;]/; it is a Mott in- 
sulator with {fia = 2,nh = 2)forU + 2zV < < 2U + 2zV, 
and so on. 

Having discussed the t = limit, now we are ready to ana- 
lyze the competition between the kinetic and potential energy 
terms of the Hamiltonian when t ^ 0. As t increases, one 
expects that the range of /i about which the ground state is 
insulating (incompressible) decreases, and that the Mott and 
CDW phases disappear at a critical value of t, beyond which 
the system becomes compressible. 

B. Transition from an Incompressible to a Compressible Phase 



A. The Atomic (t = 0) Limit 

To understand the zero-temperature (T = 0) phase diagram 
of the extended BH model given in Eq. ([T]i, we start by analyz- 
ing the atomic (t = 0) limit. In this limit, since the kinetic en- 
ergy vanishes, the boson number operator n^ commutes with 
all of the remaining terms of the Hamiltonian. Therefore, ev- 
ery lattice site is occupied by a fixed number n^ of bosons and 
the system is insulating. 

When V ~ 0, the ground-state boson occupancy is the 
same for every lattice site such that (n^) — uq where (...) 
is the thermal average, and the average boson occupancy no is 
chosen to minimize the ground-state energy for a given /i (no 
is an integer here and should not be confused with the conden- 
sate fraction of a superfluid). It turns out that the ground-state 
energy of the no state is degenerate with that of the no + 1 state 
at II — UnQ. This means that the chemical potential width of 
all Mott lobes is U, and that the boson occupancy increases 



To determine the phase boundary between the incompress- 
ible (Mott or CDW insulators) and the compressible (super- 
fluid or supersolid) phases, we need the energies of the Mott 
and CDW phases and of their defect states as a function of t. 
The defect states are characterized by exactly one extra par- 
ticle or hole which moves coherently throughout the lattice. 
At the point where the energy of the incompressible state be- 
comes degenerate with its defect state, the system becomes 
compressible assuming that the compressibility approaches 
zero continuously at the phase boundary. Therefore, the phase 
boundary between the Mott and superfluid phases is deter- 
mined by 

i^j5ottK) = i?ZttK), (2) 

E^IM = E^Uno), (3) 

where £^Mott("o) is the energy of the Mott phase with no 
bosons on every lattice site, and £^Mott('^o) and E^^^^^{no) 
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are the energies of the Mott-defect phases with exactly one 
extra particle or hole, respectively. These conditions deter- 
mine the phase boundaries of the particle and hole branches 
of the Mott insulating lobes, /i^ott ^'^'^ /^Motf respectively, as 
a function of t, U, V and uq. Similarly the phase boundary 
between the CDW and supersolid phases is determined by 



hoi 



pins 
^CDW 



CDW 



(4) 
(5) 



where E}j^t^ (n^ , n^) is the energy of the CDW phase with Ua 
and Tif, bosons on alternating lattice sites, and E^^^{na, rib) 
and i?cDW ' ""fc) '■^^ energies of the CDW-defect phases 
with exactly one extra particle or hole, respectively. These 
conditions determine the phase boundaries of the particle and 
hole branches of the CDW insulating lobes, /^cdw McdW' 
respectively, as a function of t, U, V, Ua and n^. Next, we 
calculate the energies of the Mott and CDW phases and of 
their defect states as a perturbative series in the hopping t. 



III. STRONG-COUPLING PERTURBATION THEORY 

We use the many-body version of Rayleigh-Schrodinger 
perturbation theory in the kinetic energy term lUSll to perform 
the expansion (in powers of the hopping) for the different en- 
ergies needed to carry out our analysis. The perturbation the- 
ory is performed with respect to the ground state of the system 
when the kinetic-energy term is absent. This technique was 
previously used to discuss the phase diagram of the on-site 
BH model jTt [St], and its results showed an excellent agree- 
ment with the Quantum Monte Carlo simulations (including 
the most recent numerical work ifiol I20I1 '). Here, we gener- 
alize this method to the extended BH model, hoping to de- 
velop an analytical approach which could also be as accurate 
as the numerical ones. However, we remark that our strong- 
coupling perturbation theory cannot be used to calculate the 
phase boundary between two compressible phases, e.g. the su- 
persolid to superfluid transition. In addition, we cannot even 
tell whether the compressible phase is a supersolid or a super- 
fluid. 



and B sublattices). We use, here and throughout, the index k 
to refer to all lattice sites, while the indices i and j are limited 
to the A and B sublattices, respectively. 

On the other hand, the wavefunctions of the defect states 
are determined by degenerate perturbation theory. To zeroth 
order in t, the wavefunctions for the particle-defect states can 
be written as 



I* 



par(0)x _ 
Mott / 



1 

t \ ' fMottT t 

=^ Ik Of, 



s(0)> 



k I ^ Mott / ' 



(8) 



1^ 



par(0)> 
CDW / 



fc=l 
M/2 



CDWB^tl^ins^W^^ (9) 



'^j I * CDW / ' 



where fjr is the eigenvector of the hopping matrix 
tkki with the highest eigenvalue (which is zt) such that 



Mott 



and /, 



CDWB 



is the eigenvec- 



tor of tjitij' (this matrix lives solely on the B sublat- 



Ej- j-Mott 
fc' Ifefe'/fc' 

tice) with the highest eigenvalue (which is z'^i^) such that 

highest eigenvalue of ty because the hopping matrix enters 
the Hamiltonian as —tij, and we ultimately want the lowest- 
energy states; similarly for the CDW phases, the coefficient 
of the matrix that enters the perturbation theory is neg- 
ative, so we want the highest eigenvalue again. The nor- 
malization condition requires that Efcli = 1 
Y!fll I /CDWB 1 2 ^ Y Similarly, to zeroth order in t, the 
wavefunctions for the hole-defect states can be written as 



1^ 



hol(0)> 
Mott / 



I* 



hol(0)x 
CDW/ 



^Z^Jk "fcl'l'Mott 

fc=l 
M/2 

n _ ^ — ^ 



CDWA. |,T,ins(0)> 
"il^CDW/' 



(10) 

(11) 



where fp^^^ is the eigenvector of Ej tijtji' (this ma- 
trix lives solely on the A sublattice) with the highest 
eigenvalue (which is z^t^) such that Eji' = 
z^i^/jCDWA -pjjg normalization condition requires that 



I /■CDWA|2 



1. 



A. Ground-State Wavefunctions at Zeroth Order in t 

For our purpose, we first need the ground-state wavefun- 
tions of the Mott and CDW phases and of their particle and 
hole defects when i = 0. To zeroth order in t, the Mott and 
CDW wavefunctions can be written as 

= n4^io). (6) 



fe=i " 

CD^^) = n %W%=tIO>, (7) 



I* 



where M is the number of lattice sites, and |0) is the vacuum 
state (here, we remind that the lattice is divided equally into A 



B. Ground-State Energies up to Third Order in t 

Next, we employ the many-body version of Rayleigh- 
Schrodinger perturbation theory in t with respect to the 
ground state of the system when t = 0, and calculate the en- 
ergies of the Mott and CDW phases and of their particle- and 
hole-defect states. To third order in t, the energy of the Mott 
state is obtained via nondegenerate perturbation theory and it 
is given by 



-^Mott('^o) 

M 



U 



no(no - 1) 



zV- 



fino 



no{na + 1) 



zt^ 



U-V 



Oit% 



(12) 
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which is an extensive quantity, that is £^Mott("o) is propor- 
tional to the total number of lattice sites M. The odd-order 
terms in t vanish for the d-dimensional hypercubic lattices 
considered in this manuscript. Notice that Eq. (fTzl l recovers 
the known result for the on-site BH model when V = 10,111] • 
Similarly, to third order in t, the energy of the CDW state is 
also obtained via nondegenerate perturbation theory and it can 
be written as 

^cBwi'^a, nb)_ _ ^^ na{na - 1) + nb{nb - 1) ^ ^y ^gnb 



M 
ria + rib 



Uaiub + 1) 



U [Ua - Ub-l) + V{znb - ZUa + 1) 



nb{na + 1) 



U{nb - Ua - 1) + V{zna - ZUb + 1) 



(13) 



which is also an extensive quantity, and the odd-order terms 
in t also vanish. Notice that Eq. (fTSl l reduces to Eq. (fT2b when 

Ua = rib = riQ as expected. 



The calculation of the defect state energies is more involved 
since it requires using degenerate perturbation theory. This is 
because when exactly one extra particle or hole is added to the 
Mott phase, it could go to any of the M lattice sites and all of 
those states share the same energy when t = 0. Therefore, 
for both Mott defect states with exactly one extra particle or 
hole, the initial degeneracy is of order M and it is lifted at first 
order in t. A lengthy but straightforward calculation leads to 
the energy of the Mott particle-defect state up to third order in 
t as 



■^Mott("-o) = EMottino) + Uno + zYhq - /i - (tiq + l)zt + ?io{(no + 1) 



l-z 2(1 -z) 



2z 



U 



U~2V U-V 



riQ 



2{U-V). 



>zr 



riQ 



no (no + 1){ 

4(^-2) 

U{U - 2V) ^ (f/ - V){U - 2V) 



3z + 3' 



[/2 ({/-y)2 
4(z2 -3z + 3) 



+ (no + 1) 
(no + 2 



z(l-z) 2z2-6z + 6 2z(l-z) 2{z^ - 3z + 3) 



t/2 
z- 1 



iU-V) 
z 



u{u-v) 4:{u-vy 



>zr 



iu-2vy 
+ o(t^). 



U{U-V) 



(14) 



This expression is valid for all d-dimensional hypercubic lattices, and it recovers the known result for the on-site BH model 
when y = 10, Si- To third order in t, we obtain a similar expression for the energy of the Mott hole-defect state given by 



£^Mott('^o) = E'^l,,{no) - U{no - 1) - zVna + n- n^zt + (no + 1){ 
- no(no + l)|(no + 1) 



no 



2z 



U 



2(1 -z) 



U-2V U-V 



""^^ \ze 



2{U -V). 



3z + 3 



{U-Vf 



z(l-z) 2z2-6z + 6 2z(l-z) 2(z2 - 3z + 3) 



4(z-2) 4(z^ -3z + 3) 

'U{U -2V) ^ {U -V){U -2V) 



+ {no ~ 1) 



[/2 
1 



z 



u{u-v) 4([/-y)2 



>zr 



{u~2vy 
+ o{t% 



U{U -V) 



(15) 



which also is valid for all d-dimensional hypercubic lattices, and recovers the known result for the on-site BH model when 
F = 0||7li. 

On the other hand, for d> 1 dimensions, when an extra particle or hole is added to the CDW phase, it could go to any of the 
Af/2 sites in the sublattice B or A, respectively (here, we remind that ria > rib is assumed in this manuscript). Therefore, for 
both CDW defect states with an extra particle or hole in d > 1 dimensions, the degeneracy is of order Af/2 and it is lifted at 
second order in t. This is because the states occupy one of the sublattices, and they cannot be connected by one hop, but rather 
require two hops to be connected. Another lengthy but straightforward calculation leads to the energy of the CDW particle-defect 
state up to third order in t as 



CDw('^i, Ub) = E'^l,^{na, rib) + Unb + zVria - ^J, + 



riairib + l)z 



{rig + l){nb + l)z 
U{rib ~ ria) + V{zna - zrib) 
nb{na + l)z 



Uairib + 2) 



U{na - nb-l) + V{znb - zria + 1) U {rib - na - 1) + V{zna - zrib + I) U {ria - nb - 2) + V{znb - zria + 2) 



nb{na + l)(z - 1) 



2na(nf, + l)(z-l) 



U {rib - ria - I) + V{zna - zrib) U{na - nt - 1) + V{znb - zria + 2) 



zt' 



0{t^ 



(16) 



This expression is valid for {d > 1) -dimensional hypercubic lattices. Notice that the odd-order terms in t vanish for these 
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lattices. To third order in t, we obtain a similar expression for the energy of the CDW hole-defect state given by 



JlaUbZ 



U{nb - na) + V{zna - znt) 
nb{na + l)z (na-l)(nb + l) 



U{na - rife - 1) + V{znb - zUa + 1) U{nb - ria - 1) + V{zna - zub + I) U {ria - Ub - 2) + V{znb - zUa + 2) 



nb{na + l)(z - 1) 



2na{nb + l){z - I) 



U [nb - Ha - I) + V{zna - ZUb) U{na -Ub-l) + V{znb - ZUa + 2) 

I 



zt^ + 0{t^). 



(17) 



This expression is also valid for {d > 1) -dimensional hyper- 
cubic lattices where the odd-order terms in t vanish. 

Notice that because the Mott defect states have corrections 
to first order in the hopping, while the CDW defects have cor- 
rections to second order in the hopping, the slopes of the Mott 
phase will be finite as < ^ 0, but they will vanish for the 
CDW lobes. Hence, the shape of the different types of insu- 
lating lobes are always different. 

In one dimension (d — 1), however, when exactly one ex- 
tra particle or hole is added to the CDW phase, the degener- 
acy of both of the CDW defect states is of order M and it 
is lifted at first order in t. This difference between d > 1 
and d — 1 makes one dimension unique, and it is the reason 
that the supersolid phase exists in higher dimensions but not 
in one ifisi [Till . In other words, due to this large degeneracy, 
an extra particle or hole immediately delocalizes the bosons 
in d — 1, and the crystalline order disappears. Since d — 1 
requires special attention, it will be addressed elsewhere, and 
we restrict the analysis here to higher dimensions. 

We would like to remark in passing that the energy dif- 
ference between the Mott and CDW phases with their defect 
states determine the phase boundary of the particle and hole 
branches. While all E^ZuM, E^ttM and £;Xu('^o) 
depend on the lattice size M, their difference does not. There- 
fore, the chemical potentials that determine the particle and 
hole branches, /i^ott ^^'^ /^Motf respectively, are indepen- 
dent of M at the phase boundaries. Similarly, while all 
EcDw(.^o.,nb), E^'^^{na,nb) and E^^°^^{na,nb) depend 
also on the lattice size M, their difference does not. There- 
fore, the chemical potentials that determine the particle and 
hole branches, /i^DW ^'^'^ A*cdW' respectively, are also inde- 
pendent of M at the phase boundaries. These observations 
indicate that the numerical Quantum Monte Carlo simulations 
which are based on Eqs. to (|5]l should not have too strong a 
dependence on M. It also shows that exact diagonalization on 
finite clusters of a sufficiently large size can also yield these 
expressions if properly analyzed to extract the coefficients of 
the power series. 



C. Extrapolation to Infinite order via Scaling Theory 

As a general rule, the third-order strong-coupling perturba- 
tion theory appears to be more accurate in lower dimensions. 
For this reason, an extrapolation technique to infinite order in 
t is highly desirable to determine more accurate phase dia- 



grams. Here, we propose a chemical potential extrapolation 
technique based on scaling theory to extrapolate our third- 
order power-series expansion into a functional form that is 
appropriate for the Mott and CDW lobes. 

It is known that the critical point at the tip of the Mott and 
CDW lobes has the scaling behavior of a (d + l)-dimensional 
XY model, and therefore the lobes have Kosterlitz-Thouless 
shapes for d = 1 and power-law shapes for d > 1. For the 
latter case considered in this manuscript, we propose the fol- 
lowing ansatz for the Mott and CDW lobes which includes the 
known power-law critical behavior of the tip of the lobes 



par, hoi 
/^Mott/CDW 

U 



^Mott/CDw(2;) 



±B 



Mott/CDW 



Mott/CDW 



(18) 



where AMott/CDw(2^) 
"2 - d,. 



-3 -I- and BMr.tt/o.r,w{x) 



^Mott/CDW2^ 
Mott/CDW ( 



^Mott/CDW 
CMott/CDW^;' + C^Mott/CDW^;'^ 

"Mott/CDW+Z^Mott/CDWa^+TMott/CDWa^^+'^Mott/CDWa^'^-i 

... are regular functions of x = dt/U, is the 

critical point which determines the location of the Mott and 
CDW lobes, and zv is the critical exponent for the {d + 1)- 
dimensional XY model which determines the shape of the 
Mott and CDW lobes near 2^Mott/CDW ^'i- GSI*' P^^^ 
sign corresponds to the particle branch, and the minus sign 
corresponds to the hole branch. The parameters flMott/CDW^ 
^Mott/CDW. CMott/CDW and ^Mott/CDW depend on U, V and 
no or {na,nb}, and they are determined by matching them 
with the coefficients given by our third-order expansion such 
thatylMott/CDw(2;) = (A'Mott/CDW + A'Mott/CDw)/^- To de- 
termine the U, V and uq or {ua, rib} dependence of the pa- 
rameters aMott/CDWj /^Mott/CDW^ 7Mott/CDW' (^Mott/CDW. 



X- 

of B 

,hol 



Mott/CDW ^^^^ expand the left hand side 

Dw(2;)(a:^Mott/CDW ~ ^)^'^ ~ (A'Mott/CDW ~ 

/^Mott/CDW )/^ in powers of x, and match the coefficients 
with the coefficients given by our third-order expansion. Then 
we fix zv at its well-known values such that zv w 2/3 for d — 
2 and zv — 1/2 for d > 2, and set (^Mott/CDW = to deter- 
mine ^Mott/CDW' /^Mott/CDW' 7Mott/CDW and a;^Qjj^(-,j-,^ 
self-consistently. 

Having discussed the strong coupling perturbation theory, 
next we present the ground-state phase diagrams for {d ~ 2)- 
and (d — 3)-dimensional hypercubic lattices. 
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(a) Two dimensions (V=0.1 U) 



3 
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1.8 



0.9 



ext 
s-c 



Mott{2) 



Mott(1) 



CDW(1,0) 
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x = dt/U 



0.12 



(b) Tliree dimensions (V=0.1U) 
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1.8 
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*<9=rr_ ^ 
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CDW(1,0) 
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0.04 0.08 
x = dt/U 
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FIG. 1: We show the chemical potential fj, (in units of U) versus 
X = dt/U phase diagram for (a) two- (d = 2) and (b) three- (d — 3) 
dimensional hypercubic lattices. We choose the nearest-neighbor re- 
pulsion as V = 0.1(7. The dotted lines correspond to phase bound- 
aries for the Mott insulator to superfluid and CDW insulator to su- 
persolid states as determined from the third-order strong-coupling 
perturbation theory (s-c). The circles correspond to the extrapolation 
fit (ext) discussed in the text. 



D. Numerical Results 

In Fig.[ri the results of the third-order strong-coupling per- 
turbation theory (dotted lines) are compared to those of the 
extrapolation technique (circles) when V ~ Q.IU. At t — 0, 
the chemical potential width of all Mott and CDW lobes 
are U and O.lzU, respectively where z = 2d, and that the 
ground state alternates between the CDW and Mott phases as 



a function of /i. For instance, the ground state is a vacuum 
(no = 0) for ^ < 0; it is a CDW with {ua = l,n& = 0) 
for < < O.lzU; it is a Mott insulator with (no = 1) for 
O.lzU < n< {1 + 0.1z)U; it is a CDW with {ua = 2, n^ ^ 
1) for (1 + 0.1z)U < < (1 + 0.2z)U; it is a Mott insulator 
with (no = 2) for (1 + 0.2z)U < M < (2 + 0.2z)U. 



(a) Mott critical points 



ISC W W V L ^ W NJ> W W W \ i / W MX M/ ^ y,, ^ / 

7s. /r\ /K /K 7K 7K /K A 7K 7K yf> 7F. TT. /f^ ?K — 7s 



^ 0.98 
> 

0.96 



0.94 "o^ 




0.3 0.6 
zV/U 

(b) CDW critical points 



0.12 



0.08 



3 

T3 

II 



0.04 




zV/U 

FIG. 2: (Color online) We show the critical points (location of the 
tips) Xc ~ dtc/U that are found from the chemical potential extrap- 
olation technique described in the text versus zV/U, where z — 2d. 
In Fig. (a), Xc's are scaled with their V = value; in infinite dimen- 
sions the exact critical hoppings for the Mott lobes are independent 
of V. In Fig. (b), comparing the extrapolated strong-coupling and 
exact mean-field results for the d ^ oo limit shows that the critical 
points for the CDW lobes become less accurate as V increases. This 
is because the coefficient of the O(t^) term in the power series be- 
comes very large when zV ~ 0.717, which also causes an unphysical 
decrease in Xc for zV > 0.7(7 after an initial increase. 
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As t increases from zero, the range of /i about which the 
ground state is a Mott insulator or CDW decreases, and the 
Mott insulator and CDW phases disappear at a critical value 
of t, beyond which the system becomes a superfluid near the 
Mott lobes or a supersolid near the CDW lobes. In addition, 
similar to what was found for the on-site BH model [|7|, |8t], 
the strong-coupling expansion overestimates the phase bound- 
aries, and it leads to unphysical pointed tips for all Mott and 
CDW lobes. This is not surprising since a finite-order per- 
turbation theory cannot describe the physics of the tricritical 
point correctly. 

In Fig. 12] we show the critical points (location of the tips) 
Xc — dtc/U versus zV/U. In Fig. |2la), Xc's of the Mott 
lobes are scaled with their V" = value. The critical points 
are calculated with the chemical potential extrapolation tech- 
nique that is based on the scaling theory with the exponent zv 
fixed to its known value. It is expected that the locations of the 
tips of the CDW lobes to increase as a function of V, because 
the presence of a nonzero V is what allowed these states to 
form in the first place (the Mott insulator critical points tend to 
move in as V increases). Comparing the extrapolated strong- 
coupling and exact mean-field (to be discussed below) results 
for the d ^ oo limit shows that the critical points for the CDW 
lobes become less accurate as V increases. It turns out that 
the coefficient of the 0{t'^) term in the power series is gener- 
ally small for the Mott lobes, but it can become very large for 
the CDW lobes when zV ^ U. We remind that we assume 
U > zV in this manuscript. As shown in Fig.|2b), This also 
causes an unphysical decrease in Xc for zV > 0.7U after an 
initial increase. Therefore, inclusion of the 0(t*) terms in the 
expansion are necessary to improve the accuracy of the phase 
boundaries near the tips of the CDW lobes when zV ~ U. 
In addition, we present a short list of V/U versus the critical 
points Xc = dtc/U in Table U for {d — 2)- and {d — 3)- 
dimensional lattices. 

As a further check of the accuracy of our perturbative ex- 
pansion, next we compare d ^ oo limit of our results to the 
mean-field one which corresponds to the exact solution on an 
{d — ^ cx3)-dimensional hypercubic lattice. 



IV. MEAN-FIELD DECOUPLING THEORY 

In the large-dimensional case, mean-field theory becomes 
exact, so examining the mean-field theory for the extended BH 
model provides another way to validate the strong-coupling 
expansion and to test to see how well the scaling result pro- 
duces the correct phase diagram. 

In constructing the mean-field theory, one first defines 
the superfluid order parameter as (pk = (bk) where (...) is 
the thermal average, and then replaces the operator bk with 
'fik + 5bk in the hopping term of Eq. ([T]i. This approximation 
decouples the two-particle hopping term into single-particle 
ones, and the resultant mean-field Hamiltonian can be solved 
via exact diagonalization in a power series of ^pk- The order 
parameter is finite {ifk 7^ 0) for the superfluid and supersolid 



ground states, and it vanishes (ipk — 0) for the Mott and CDW 
phases. Therefore, (fk 0+ signals the phase boundary be- 
tween an incompressible and a compressible phase. The gen- 
eralized order parameter equation to the case of 7^ can be 
written as 1I21I1 



nk 



1 



nk 



Uuk + P - M U{nk -l) + V, 



dip 



(19) 

where ipk = J2{k') 'fik' is the sum of the order parameters at 
sites k' neighboring to site k, and V^'^''' — V J2{ki) ^^k' is the 
interaction of one atom with sites k' neighboring to the site k. 

To determine the phase boundary between the Mott and su- 
perfluid phases from Eq. ( fT9] l, we set tp^ — ifo, (pk = z(po, 



and K 



dip 



zVuq. Since ipo 0+ near the phase boundary. 



Eq. (figj l can be satisfied only if 



no + 1 



no 



1 _ 

zt IIuq + zVhq — /i U {no — 1) + zVno — /i 



, (20) 



which gives a quadratic equation for pi. Notice that this equa- 
tion recovers the known result for the on-site BH model when 
V = [l6l|22t], and it can be easily solved to obtain 



par, hoi 
A^Mott 



U{no - l/2) + zVno - zt/2 
± ^/U^/4: - U{no + l/2)zt + zH"^, 



(21) 



where the plus sign corresponds to the particle branch, and 
the minus sign corresponds to the hole branch. In the d — > 00 
limit, we checked that our strong-coupling perturbation results 
for the Mott lobes agree with this exact solution when the lat- 
ter is expanded out to third order in t, providing an indepen- 
dent check of the algebra (one must note that the terms V and 
2V that appear in the denominator vanish in the limit when 
d ^ 00 because V oc 1/d). Equation ( 1211 1 also shows that the 
Mott lobes are separated by zV, but their shapes are indepen- 
dent of V; in particular, the critical points for the Mott lobes 
are independent of V. 

To determine the phase boundary between the CDW and 
supersolid phases from Eq. ( fT9] l, we set (pi — cpA, 'fi — zipB 
and V^^^ = zVnB for i G A sublattice, and we set ipj — (ps, 
Cpj — zip A and V^'^''' — zVnA for j G B sublattice. This leads 
to two coupled equations for pA and pB- Since {p^A, fs} 
0+ near the phase boundary, Eq. (fT9] l can be satisfied only if 



1 



1 



Una + zVnb 
nh + 1 



U{na - 1) + zVnb - pL 
rib 



Uub + zVna — pi U {nb — 1) + zVna — fi 



(22) 



which gives a quartic equation for pi. Since a simple closed 
form analytic solution for pi is not possible, we solve Eq. (l22l l 
with Mathematica for each of the CDW lobes separately. In 
the d 00 limit, we also checked that our strong-coupling 
perturbation results for the CDW lobes agree with this exact 
solution when the latter is expanded out to third order in t, 
providing again an independent check of the algebra. 
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TABLE I: We list the critical points (location of the tips) Xc = dtc/U that are found from the chemical potential extrapolation technique 
described in the text. 



v/u 


Two dimensions 


Three dimensions 


CDW(1,0) 


Mott(l) 


CDW(2,1) Mott(2) 


CDW(1,0) 


Mott(l) 


CDW(2,1) Mott(2) 


0.00 


- 


0.117 




0.0691 


- 


0.0981 




0.0578 


0.01 


0.00929 


0.117 


0.00465 


0.0689 


0.0143 


0.0977 


0.00717 


0.0576 


0.02 


0.0183 


0.116 


0.00916 


0.0687 


0.0278 


0.0974 


0.0139 


0.0574 


0.03 


0.0270 


0.116 


0.0135 


0.0684 


0.0405 


0.0970 


0.0203 


0.0571 


0.04 


0.0354 


0.116 


0.0178 


0.0682 


0.0522 


0.0966 


0.0263 


0.0569 


0.05 


0.0434 


0.115 


0.0219 


0.0680 


0.0630 


0.0962 


0.0317 


0.0567 


0.06 


0.0512 


0.115 


0.0258 


0.0678 


0.0723 


0.0958 


0.0367 


0.0564 


0.07 


0.0586 


0.115 


0.0295 


0.0676 


0.0814 


0.0955 


0.0411 


0.0562 


0.08 


0.0656 


0.114 


0.0331 


0.0673 


0.0888 


0.0951 


0.0449 


0.0559 


0.09 


0.0721 


0.114 


0.0365 


0.0671 


0.0947 


0.0947 


0.0480 


0.0557 


0.10 


0.0783 


0.114 


0.0396 


0.0669 


0.0990 


0.0942 


0.0502 


0.0555 



oo dimensions (zV=0.4U) 





m-f 




ext • 




s-c 


iviott(2) J 




CDW(2J^-:-. 




Mott(1) 


CDW(1,0) 



0.04 0.08 0.12 

x = dt/U 

FIG. 3: (Color online) We show the chemical potential /i (in units of 
U) versus x — dt/U phase diagram for a (d ^ oo)-dimensional hy- 
percubic lattice. Here the nearest-neighbor repulsion scales inversely 
with d such that zV — 0.4(7. The dotted lines correspond to phase 
boundaries for the Mott insulator to superfluid and CDW insulator to 
supersolid states as determined from the third-order strong-coupling 
perturbation theory (s-c). The circles correspond to the extrapolation 
fit (ext) discussed in the text. The red solid lines correspond to phase 
boundaries for the Mott insulator to superfluid and CDW insulator 
to supersolid states as determined from the mean-field theory (m-f) 
which becomes exact for d — + cxj. 



In Fig. |3] the results of the third-order strong-coupling per- 
turbation theory (dotted lines) is compared to those of the ex- 
act mean-field theory (red solid lines) and of the extrapolation 
technique (circles) for an infinite (d cxD)-dimensionaI hy- 
percubic lattice when zV = OAU. Notice that, in infinite 
dimensions, both t and V must scale inversely with d such 



that dt and dV are finite. The extrapolated solutions are indis- 
tinguishable from the exact ones for the Mott lobes, and they 
are within 5% of each other for the tips of the CDW lobes. It 
turns out that this minor disagreement around the tips of the 
CDW lobes is due to the large coefficient of the 0{t^) term 
in the power-series expansion. Therefore, we conclude that, 
even in infinite dimensions, the agreement of the third-order 
strong-coupling perturbation theory with the exact mean-field 
theory is quite good. 

V. CONCLUSIONS 

We analyzed the zero temperature phase diagram of the ex- 
tended Bose-Hubbard (BH) model with on-site and nearest- 
neighbor boson-boson repulsions in (d > l)-dimensionaI 
hypercubic lattices. We used the many-body version of 
Rayleigh-Schrodinger perturbation theory in the kinetic en- 
ergy term with respect to the ground state of the system when 
the kinetic energy term is absent. This technique was pre- 
viously used to discuss the phase diagram of the on-site BH 
model \2b iSJ, and its extrapolated results showed an excel- 
lent agreement with the recent Quantum Monte Carlo simu- 
lations lfT9ll20|] . Here, we generalized this method to the ex- 
tended BH model, hoping to develop an analytical approach 
which could be as accurate as the numerical ones. 

We derived analytical expressions for the phase bound- 
aries between the incompressible (Mott or charge-density- 
wave (CDW) insulators) and compressible (superfluid or su- 
persolid) phases up to third order in the hopping t. However, 
we remark that the strong-coupling perturbation theory devel- 
oped here cannot be used to calculate the phase boundary be- 
tween two compressible phases, e.g. the supersolid to super- 
fluid transition. We also proposed a chemical potential extrap- 
olation technique based on the scaling theory to extrapolate 
our third-order power series expansion into a functional form 
that is appropriate for the Mott or CDW lobes. 

We believe some of our results could potentially be ob- 
served with ultracold dipolar Bose gases loaded into optical 
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lattices ll2ll I23I1 . This is motivated by the recent success in 
observing superfluid to Mott insulator transition with ultra- 
cold point-like Bose gases loaded into optical lattices. Such 
lattices are created by the intersection of laser fields, and 
they are nondissipative periodic potential energy surfaces for 
the atoms. An ultracold dipolar Bose gas can be realized in 
many ways with optical lattices. For instance, heteronuclear 
molecules which have permanent electric dipole moments, 
Rydberg atoms which have very large induced electric dipole 
moment, or Chromium-like atoms which have large intrinsic 
magnetic moment, etc. can be used to generate sufficiently 
strong long-ranged dipole-dipole interactions. 

This work can be extended in several ways if desired. For 
instance, our current results for the CDW phase are not di- 
rectly applicable to the one-dimensional case. We are cur- 
rently working on this problem and will report our results 
elsewhere. In addition, it turns out that the coefficient of the 
O(t^) term in the power series is generally small for the Mott 
lobes, but it can become very large for the CDW lobes when 



zV ^ U. Therefore, inclusion of the O(t^) is necessary to im- 
prove the accuracy of the phase boundaries near the tips of the 
CDW lobes when zV ~ U. Lastly, one can include the next- 
nearest-neighbor repulsion term to the current model, which 
would lead to additional CDW phases. One can also examine 
how the momentum distribution changes with the hopping in 
the CDW phase, or in the Mott phase when there is a nearest- 
neighbor repulsion. This last calculation could have direct rel- 
evance for experiments on these systems and would generalize 
recent results for theV = case L24i] . 
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